function [valueCI,error]=bootstrap(totalData,value,err,E0,Kb,R,h,webpage)

tic
%setting some needed constants
nbs=1000;%number of samples to use (should be somewhere around 500-1000)
valueCI=zeros(nbs,5);

%all the available data
            [n,m]=size(totalData);
            Temp = totalData(1,2:m);                     %All temperature values
            Time = totalData(2:n,1);                     %All time values        
            [TT,tt]=meshgrid(Temp,Time);            %Make a grid of all possibilities
            TT=reshape(TT,(n-1)*(m-1),1);           %reshape grid into vector
            tt=reshape(tt,(n-1)*(m-1),1);  
            yy=reshape(totalData(2:n,2:m),(n-1)*(m-1),1);%all data values in one vector
            nyy=length(yy);
            
%work out nbs subsamples, each equal to 1/10th of the whole dataset
n=round(nyy/10);
i=ceil(rand(n*nbs,1).*nyy);  %randomnly picks the data values that will be included.
    clear y t T
    y(1:n,1:nbs)=reshape(yy(i),n,nbs);%creates a matrix with all the random subsets (one per column)
    t(1:n,1:nbs)=reshape(tt(i),n,nbs);
    T(1:n,1:nbs)=reshape(TT(i),n,nbs);


    
    waith = waitbar(0,'Please wait');
for q=1:nbs %work through each of the subsets finding the best-fit
    ss=toc./60
    waitbar(q/nbs,waith,['Process ',num2str((q/nbs*100)),' % complete in ',num2str(ss), ' minutes'])
%reset the values back to the optimum
valueCI(q,1:5)=value(1,:);

for jj=1:5  %you can change this range to make a better fit which is also a little slower.
%make up a 4 dimensional value grid
nx=2;   %making this number larger will make the fit alot slower.
step=0.1.*10^(-jj/2);%creating a generic step size
r=step.*valueCI(q,:);%creating a range to look for each variable
r(3)=r(3)*2;r(1)=r(1)/2;r(4)=r(4)/2;%adjusting the range a bit for each variable
dx=r./nx;%creating a step for each variable
xx1=valueCI(q,1)-r(1):dx(1):valueCI(q,1)+r(1);
xx2=valueCI(q,2)-r(2):dx(2):valueCI(q,2)+r(2);
xx3=valueCI(q,3)-r(3):dx(3):valueCI(q,3)+r(3);
xx4=valueCI(q,4)-r(4):dx(4):valueCI(q,4)+r(4);
[x1,x2,x3,x4]=ndgrid(xx1,xx2,xx3,xx4);
[nn,mm,oo,pp]=size(x1);

%create vectors of all possible values to try
xx1=reshape(x1,nn*mm*oo*pp,1);           
xx2=reshape(x2,nn*mm*oo*pp,1);           
xx3=reshape(x3,nn*mm*oo*pp,1);           
xx4=reshape(x4,nn*mm*oo*pp,1);           

%work out the error for each of these possibilities
   ERR=mean((calc(T(:,q),t(:,q),xx1,xx2,xx3,xx4,E0,h,Kb,R)-y(:,q)*ones(1,nn*mm*oo*pp)).^2)';   
%find the best one
   j=find(ERR==min(ERR));
%store the best one for the next search at a finer resolution or step size
   valueCI(q,1)=xx1(j);
   valueCI(q,2)=xx2(j);
   valueCI(q,3)=xx3(j);
   valueCI(q,4)=xx4(j);
   valueCI(q,5)=ERR(j);
   %plot(ERR)
   %pause
end   
end

%every sub-sample has been processed, work out the error characteristics.
error=zeros(4,5);
if nbs>=5
    error(1,:)=std(valueCI);    %standard deviation
    error(2,:)=mean(valueCI);   %mean
    n5=round(nbs.*0.05);       
    n95=round(nbs.*0.95);
    v=sort(valueCI);            %Sorting all the values
    error(3,:)=v(n5,:);         %fifth percentile
    error(4,:)=v(n95,:);        %ninety-fifth percentile
    
if webpage
    figure(2)
    clf
    subplot(221)
        hist(valueCI(:,1))
        title(['Variability around EACAT'])
        hold on
        plot([value(1,1) value(1,1)],[0 nbs/3],'r')
        plot([error(2,1) error(2,1)],[0 nbs/3],'k')
        plot([error(3,1) error(3,1)],[0 nbs/3],'k--')
        plot([error(4,1) error(4,1)],[0 nbs/3],'k--')
        legend('sub-samples','original fit','new fit','confidence interval')
        ylabel(['Frequency'])
    subplot(222)
        hist(valueCI(:,2))
        title(['Variability around EAl'])
        hold on
        plot([value(1,2) value(1,2)],[0 nbs/3],'r')
        plot([error(2,2) error(2,2)],[0 nbs/3],'k')
        plot([error(3,2) error(3,2)],[0 nbs/3],'k--')
        plot([error(4,2) error(4,2)],[0 nbs/3],'k--')
        legend('sub-samples','original fit','new fit','confidence interval')
        ylabel(['Frequency'])
    subplot(223)
        hist(valueCI(:,3))
        title(['Variability around HEQ'])
        hold on
        plot([value(1,3) value(1,3)],[0 nbs/3],'r')
        plot([error(2,3) error(2,3)],[0 nbs/3],'k')
        plot([error(3,3) error(3,3)],[0 nbs/3],'k--')
        plot([error(4,3) error(4,3)],[0 nbs/3],'k--')
        legend('sub-samples','original fit','new fit','confidence interval')
        ylabel(['Frequency'])
    subplot(224)
        hist(valueCI(:,4))
        title(['Variability around TEQ'])
        hold on
        plot([value(1,4) value(1,4)],[0 nbs/3],'r')
        plot([error(2,4) error(2,4)],[0 nbs/3],'k')
        plot([error(3,4) error(3,4)],[0 nbs/3],'k--')
        plot([error(4,4) error(4,4)],[0 nbs/3],'k--')
        legend('sub-samples','original fit','new fit','confidence interval')
        ylabel(['Frequency'])
     figure(3)
     clf
     plot(valueCI(:,5),'.')
     ERR=mean((calc(TT,tt,value(1,1),value(1,2),value(1,3),value(1,4),E0,h,Kb,R)-yy).^2)';
     hold on
     plot([0 nbs],[ERR ERR],'r')
     ERR=mean((calc(TT,tt,error(2,1),error(2,2),error(2,3),error(2,4),E0,h,Kb,R)-yy).^2)';
     plot([0 nbs],[ERR ERR],'k--')
     ylabel(['Residual mean square error'])
     legend('subsamples','original fit','new fit')
     title(['Error analysis'])
end
end
close(waith)
toc
return

   